function C = windowed_est(x, p)
	[N,n] = size(x);
	H = [];
	for i=1:n
		H = [H; transpose(toeplitz([x(1,i); zeros(p,1)], [x(:,i); zeros(p,1)]))];
	end
	H = reshape(H(:), N+p, [])';
	C = 1/N * H * H';
	C = (C+C')/2; % avoid asymmetry due to floating point error
end
